 classdef etsEstim
    %etsEstim Estimates of model for pollution abatement
    %   
    
    properties
        
        % Data objects used in the model
        plants;   % Plant-level data on characteristics
        bids;     % Plant X compliance period X bid level data on bids offered
        N;
        
        % Simulation parameters and bids selected
        S = 200;  % Number of simulations 
        
        % Estimation options
        modnum;
        models = {'Random effects','Plant FEs','Plant-period FEs','Plant-period FEs + Hetero', 'Plain Plant-period FEs' };
        sample = 'First half';
        specs;
        fits;
        
        % Variable names for regression specifications
        fe_periods;
        fe_plants;
        fe_pp;
        
        % Parameters of emissions rate distribution to be fit in control
        parER;
        
        % Initialize estimates of abatement parameters
        alpha;  % Pollution level and elasticity wrt heat output
        beta;   % Elasticity of MAC wrt emissions
        
        % Fits to regulation in control group
        cmodels = {'Constant','Heat output'};
        cspecs  = {'logERit ~ 1', 'logERit ~ logHi'};
        cfit;
        
    end
    
    methods
        
        function obj = etsEstim( plant_data, bid_data )
            % Initialize etsEstim object
            
            % Store plant and bidding data
            obj.plants = plant_data;
            obj.bids   = bid_data; 
            obj.N      = size(obj.bids,1);
            
        end
        
        function obj = estimate( obj )
            % Estimate parameters of ETS model
            
            switch obj.models{obj.modnum}
                
                case 'Random effects'
                    obj = obj.estimateRandomEffects;
                    
                case 'Plant FEs'
                    obj = obj.estimatePlantEffects;
                    
                case 'Plant-period FEs'
                    obj = obj.estimatePlantPeriodEffects;

                case 'Plant-period FEs + Hetero'
                    obj = obj.estimatePlantPeriodHeteroEffects;
                
                case 'Plain Plant-period FEs'
                    obj = obj.estimatePlainPlantPeriodEffects;
            end
                
        end
                
        function obj = preEstimation( obj )
            
            % Transform data as needed for estimation
            obj.bids.log_bid_p = log(obj.bids.bid_p);
            obj.bids.logH      = log(obj.bids.Hi);
            obj.bids.logE      = log(obj.bids.bid_qE);
        
            % Mark sample as bids with non-missing price and quantity
            require = {'log_bid_p','logH','logE'};
            obj.bids.insample = ...
                ~any( ismissing(obj.bids(:,require),{NaN Inf -Inf}), 2);
            
            switch obj.sample
                case 'First half'
                    obj.bids.insample = obj.bids.insample & ...
                        obj.bids.D_first_half == 1;
            end
            
            % Add period dummies
            if isempty( obj.fe_periods )
                [ obj.bids, obj.fe_periods ] = ...
                    dummyOut( obj.bids, 'id_period', 1 );
            end
                        
            switch obj.models{obj.modnum}
                
                case 'Random effects'
                                
                case 'Plant FEs'
                    [ obj.bids, obj.fe_plants ] = ...
                        dummyOut( obj.bids(obj.bids.insample,:), 'id_plant', [] );
                    
                case 'Plant-period FEs'
                    obj.bids.id_pp = obj.bids.id_period * 10000 + ...
                        obj.bids.id_plant;
                    obj.bids = movevars(obj.bids,'id_pp','After','id_period');
                    [ obj.bids, obj.fe_pp ] = ...
                        dummyOut( obj.bids(obj.bids.insample,:), 'id_pp', [] );

                case 'Plant-period FEs + Hetero'
                    if sum(strcmp('id_pp',obj.bids.Properties.VariableNames)) == 0
                        obj.bids.id_pp = obj.bids.id_period * 10000 + ...
                            obj.bids.id_plant;
                        obj.bids = movevars(obj.bids,'id_pp','After','id_period');
                        [ obj.bids, obj.fe_pp ] = ...
                            dummyOut( obj.bids(obj.bids.insample,:), 'id_pp', [] );
                    end
                    % Construct interaction between logE and apcd_max indicators
                    obj.bids.logE_1 = obj.bids.logE .* (obj.bids.apcd_max == 1); % logE x 1[apcd_max = "bag filter"]
                    obj.bids.logE_2 = obj.bids.logE .* (obj.bids.apcd_max == 2); % logE x 1[apcd_max = "cyclone"]
                    obj.bids.logE_12 = obj.bids.logE .* ((obj.bids.apcd_max == 1) | (obj.bids.apcd_max == 2)); % logE x 1[apcd_max = "cyclone"/"bag filter"]
                    obj.bids.logE_3 = obj.bids.logE .* (obj.bids.apcd_max == 3); % logE x 1[apcd_max = "scrubber"]
                    obj.bids.logE_4 = obj.bids.logE .* (obj.bids.apcd_max == 4); % logE x 1[apcd_max = "ESP"]
                    obj.bids.logE_34 = obj.bids.logE .* ((obj.bids.apcd_max == 3) | (obj.bids.apcd_max == 4)); % logE x 1[apcd_max = "scrubber"/"ESP"]
                
                    % Not needed as dummies already created in earlier run
                % case 'Plain Plant-period FEs'
                %     obj.bids.id_pp = obj.bids.id_period * 10000 + ...
                %         obj.bids.id_plant;
                %     obj.bids = movevars(obj.bids,'id_pp','After','id_period');
                %     [ obj.bids, obj.fe_pp ] = ...
                %         dummyOut( obj.bids(obj.bids.insample,:), 'id_pp', [] );
            end
            
        end
         
        function obj = estimateRandomEffects( obj )
            % Estimate parameters of ETS model with plant random effects
           
            % Add necessary variables to the dataframe
            obj = obj.preEstimation;
            
            % Matlab example of fitlme specification
            %   'CityMPG ~ Horsepower + (1|EngineType) + (Horsepower-1|EngineType)'
            %           'Fixed effects'            'Random effects'
            % In matlab terminology, an explanatory variable with a
            %   constant coefficient has a fixed effect.
            
            % Build specification
            depvar    = 'log_bid_p ~ ';
            periodFEs = strjoin(obj.fe_periods,' + ');
            X1 = ['logH + logE + ' periodFEs ]; % Variables with constant effects
            X2 = ['+ (1|id_plant)'];            % Variables with random effects
            obj.specs{1} = [ depvar X1 X2 ];
       
            % Fit specification
            obj.fits{1} = fitlme(obj.bids(obj.bids.insample,:),obj.specs{1});
        end        
        
        function obj = estimatePlantEffects( obj )
            % Estimate parameters of ETS model with plant fixed effects
           
            % Add necessary variables to the dataframe
            obj = obj.preEstimation;
             
            % Build specification
            depvar    = 'log_bid_p ~ ';
            periodFEs = strjoin(obj.fe_periods,' + ');
            plantFEs  = strjoin(obj.fe_plants,' + ');
            X1 = ['-1 + logE + ' periodFEs ' + ' plantFEs ]; % Constant effects
            X2 = [''];                                       % Random effects
            obj.specs{2} = [ depvar X1 X2 ];
       
            obj.fits{2} = fitlme(obj.bids(obj.bids.insample,:),obj.specs{2});
        end        
        
        function obj = estimatePlantPeriodEffects( obj )
            % Estimate parameters of ETS model with plant-period fixed effects
           
            % Add necessary variables to the dataframe
            obj = obj.preEstimation;
            
            % Build specification
            depvar = 'log_bid_p ~ ';
            ppFEs  = strjoin(obj.fe_pp,' + ');
            X1 = ['-1 + logE + ' ppFEs ]; % Constant effects
            X2 = [''];                    % Random effects
            obj.specs{3} = [ depvar X1 X2 ];
       
            obj.fits{3} = fitlme(obj.bids(obj.bids.insample,:),obj.specs{3});
        end    

        function obj = estimatePlantPeriodHeteroEffects( obj )
            % Estimate parameters of ETS model with plant-period fixed effects 
            % + heterogeneity by abatement technology
           
            % Add necessary variables to the dataframe
            obj = obj.preEstimation;
            
            % Build specification
            depvar = 'log_bid_p ~ ';
            ppFEs  = strjoin(obj.fe_pp,' + ');
            % X1 = ['-1 + logE_1 + logE_2 + logE_3 + logE_4 + ' ppFEs ]; % Constant effects
            % X1 = ['-1 + logE_1 + logE_2 + logE_34 + ' ppFEs ]; % Constant effects
            X1 = ['-1 + logE_12 + logE_34 + ' ppFEs ]; % Constant effects
            X2 = [''];                    % Random effects
            obj.specs{4} = [ depvar X1 X2 ];
            
            obj.fits{4} = fitlme(obj.bids(obj.bids.insample,:),obj.specs{4});
        end    
        
        function obj = estimatePlainPlantPeriodEffects( obj )
            % Estimate parameters of ETS model with plant-period fixed effects
           
            % Add necessary variables to the dataframe
            obj = obj.preEstimation;
            
            % Build specification
            depvar = 'bid_p ~ ';
            ppFEs  = strjoin(obj.fe_pp,' + ');
            X1 = ['-1 + ' ppFEs ]; % Constant effects
            X2 = [''];                    % Random effects
            obj.specs{5} = [ depvar X1 X2 ];
       
            obj.fits{5} = fitlme(obj.bids(obj.bids.insample,:),obj.specs{5});
        end  

        %% Fit emissions rate in control group
        function obj = estimateControlStringency( obj )
            % Estimate command and control stringency in control group
            rng(0);
            
            % Sample restriction
            keep = ~isnan(obj.plants.ERit) & obj.plants.treatment==0;
            T = max(obj.plants.id_period);
            obj.parER = table('Size',[T 2],...
                'VariableNames',{'mu' 'sigma'},...
                'VariableTypes',{'double','double'});
            
            % Initialize model output
            obj.cfit = cell(T,2);
            
            % Create variables for regression of emissions rates
            obj.plants.logERit = log(obj.plants.ERit);
            obj.plants.logHi   = log(obj.plants.Hi);
            
            for t = 1:T
                
                % Distribution of emissions rates
                theta = lognfit(obj.plants.ERit(keep & obj.plants.id_period == t));
                obj.parER.mu(t) = theta(1);
                obj.parER.sigma(t) = theta(2);
                
                % Regression of emissions rate on constant
                obj.cfit{t,1} = fitlme(...
                    obj.plants(keep & obj.plants.id_period == t,:),...
                    obj.cspecs{1});
                
                % Regression of emissions rate on heat output capacity
                obj.cfit{t,2} = fitlme(...
                    obj.plants(keep & obj.plants.id_period == t,:),...
                    obj.cspecs{2});
            end
            
            fprintf(1,'Estimated emissions rate distribution parameters\n');
            display(obj.parER);
        end

        function obj = etsAbatementModelsWrapper( obj )
            rng(0);
            
            %% Estimate ETS abatement model in treatment
            % Model 1: plant random effects
            obj.modnum = 1;
            obj = obj.estimate;
            display(obj.fits{1});

            % Model 2: plant fixed effects
            obj.modnum = 2;
            obj = obj.estimate;
            display(obj.fits{2});

            % Model 3: plant-period fixed effects
            obj.modnum = 3;
            obj = obj.estimate;
            display(obj.fits{3});

            % Plant-period fixed effects + Hetero
            obj.modnum = 4;
            obj = obj.estimate;
            display(obj.fits{4});

            % Model 5: plain plant-period fixed effects (ie step function)
            obj.modnum = 5;
            obj = obj.estimate;
            display(obj.fits{5});
            
        end
    
    end % End methods block
    
 end
 
 function [ table, names ] = dummyOut( table, dcateg, omit )
     % Dummy out variable dcateg in data table table
     %   Create variables dcateg_n for each category of dcateg
     if nargin < 3
         omit = [];
     end
     
     % Record unique values of dcateg variable
     [ c ] = unique( table{:,{dcateg}} );
     names = {};
     
     % Omit a value if desired
     if ~isempty(omit)
         c = c( ~ismember(c,omit) );
     end
     
     for i = length(c):-1:1
         dcatname = [ dcateg '_' num2str(c(i)) ];
         table = addvars(table, int8(table{:,{dcateg}} == c(i)) ,...
             'After',dcateg,'NewVariableNames',dcatname);
         names{i} = dcatname;
     end
     
 end
 
 
 
